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FINITE SIZE CORRECTIONS IN RANDOM MATRIX THEORY 
AND ODLYZKO’S DATA SET FOR THE RIEMANN ZEROS 

PETER J. FORRESTER AND ANTHONY MAYS 


Abstract. Odlyzko has computed a data set listing more than 10® successive 
Riemann zeros, starting at a zero number beyond 10®^. The data set relates 
to random matrix theory since, according to the Montgomery-Odlyzko law, 
the statistical properties of the large Riemann zeros agree with the statistical 
properties of the eigenvalues of large random Hermitian matrices. Moreover, 
Keating and Snaith, and then Bogomolny and collaborators, have used N x N 
random unitary matrices to analyse deviations from this law. We contribute 
to this line of study in two ways. First, we point out that a natural process to 
apply to the data set is to thin it by deleting each member independently with 
some specified probability, and we proceed to compute empirical two-point 
correlation functions and nearest neighbour spacings in this setting. Second, 
we show how to characterise the order l/N'^ correction term to the spacing 
distribution for random unitary matrices in terms of a second order differen¬ 
tial equation with coefficients that are Painleve transcendents, and where the 
thinning parameter appears only in the boundary condition. This equation 
can be solved numerically using a power series method. Comparison with the 
Riemann zero data shows accurate agreement. 


1. Introduction 

The celebrated Riemann hypothesis asserts that all the non-trivial zeros of the 
Riemann zeta function C(s) are on the so-called critical line, Res = 1/2. We will 
refer to zeros of C(s) on the critical line as Riemann zeros, and restrict attention 
to those in the upper half plane — those in the lower half plane are obtained 
by complex conjugation. There is an intriguing relationship between the Riemann 
zeros and random matrix theory. It originates from the viewpoint that the Riemann 
zeros, for distances far up the critical line, are best considered in a statistical sense 
rather than as a deterministic sequence. The first step in a statistical analysis is to 
scale the sequence so that locally the density of the zeros far up the critical line is 
unity. This is straightforward from knowledge of the fact that at position ^ iE 
along the critical line, and with E :$> the density {p say) is given by [55l pg. 280] 

Montgomery |U10^ made a study of the density of the scaled zeros about a fixed 
zero, with the latter averaged over some window of size AE, with AE <C E. This 
quantity, referred to as the two-point correlation function, is probed by computing 
the average value of a test function /{E — Eq), with Eq corresponding to the fixed 
zero. Subject to the technical assumption that the class of test function used has 
Fourier transform supported on |fc| < 2tt, the limiting two-point correlation function 
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was proved to be equal to 


lim pf {E,E + s) = l-r-^)\ 

Af—>oo ' ' \ TTS / 


( 1 . 2 ) 


where the superscript ‘R’ denotes that this quantity pertains to the Riemann zeros 
data. On the other hand, let X be an iV x iV random matrix with standard complex 
Gaussian entries, and form the Hermitian matrix ^(X + Xl) — the set of such 
matrices specifies the Gaussian Unitary Ensemble (GUE). After normalisation so 
that the eigenvalues in the bulk of the spectrum have unit density, the N ^ oo form 


of the two-point correlation function is precisely (1.2); see e.g. m Prop. 7.1.1]. 

This coincidence, first observed by Dyson (see e.g. [TH pg. 159]), is in keeping 
with and extends the so-called Hilbert-Poly a conjecture |56j that the Riemann zeros 
correspond to the eigenvalues of some unbounded self-adjoint operator. Gonsider- 
ations in quantum chaos (see e.g. mm give evidence that generically the large 
eigenvalues of such operators have local statistics coinciding with the eigenvalues of 
large random Hermitian matrices. The latter must be constrained to have real en¬ 
tries should there be a time reversal symmetry. Thus this line of reasoning suggests 
that asymptotically the Riemann zeros have the same local statistical properties as 
the large eigenvalues of a chaotic quantum system without time reversal symmetry. 
This assertion is essentially a statement of what now is called the Montgomery- 
Odlyzko law or the GUE conjecture; see also [111] the reviews [IlllS], and the 
thesis [KH] . 

At an analytic level, further evidence for the validity of the law was given by 
Rudnick and Sarnak [51] , who extended Montgomery’s result to the general fc-point 
correlation function by obtaining the explicit functional form 


det 


sin7r(a;j — xi) 

TT{Xj - Xl) 


(1.3) 




again subject to a constraint on the Fourier transform of the class of test functions 
involved; see the recent work [52| for some refinements. By use of a conjecture of 
Hardy and Littlewood relating to the pair correlation of prime numbers, and also 
at the expense of introducing some non-rigorous working, Bogomolny and Keating 


piiniiii] removed this constraint. The /c-point correlation ( |1.3[ ) is precisely that 
for the eigenvalues of large Hermitian random matrices in the bulk scaling limit; 
see again Prop. 7.1.1]. At a numerical level, Odlyzko has made high precision 
computations of the 10^®-th Riemann zero and over 70 million of its neighbours 
[46] . and (later) of the 10^^-nd zero and one billion of its neighbours [47]. (The 
zeros on the critical line are conventionally numbered in the order of their distance 
from the real axis, with the first zero at Ims « 14.134725 |S31 A058303].) This 
data exhibits consistency with random matrix theory for statistical quantities such 
as the pair correlation function outside the range known rigorously from the work 
of Montgomery, the distribution of the spacing between neighbouring zeros, the 
variance of the fluctuation of the number density in an interval etc.; for a popular 
account of this line of research see |35) . 

Notwithstanding the extraordinary distance along the critical line achieved in 
Odlyzko’s calculations, it turns out that finite size corrections to the limiting be¬ 


haviour occur on the scale of the logarithm of the distance as suggested by (1.1), 


and thus are of significance in the interpretation of the data from the viewpoint of 
the random matrix predictions. In addressing the question of the functional form of 
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the finite size corrections from the viewpoint of random matrix theory, Keating and 
Snaith [33] (see also the review [30] and the later work m) were led to a totally 
unexpected conclusion: the correct model for this purpose is not complex random 
Hermitian matrices, but rather N x N random unitary matrices chosen with Haar 


measure, and N related to E so as to be consistent with (1.1). Random unitary 
matrices with Haar measure result, for example, by applying the Gram-Schmidt 
orthonormalisation procedure to a matrix of standard complex Gaussians. The 
statistical quantities considered in [33] relate to the value distribution of the loga¬ 
rithm of the Riemann zeta function on the critical line, which for random matrices 
corresponds to the value distribution of the characteristic polynomial. Coram and 
Diaconis [23] subsequently used the random unitary matrix model to predict the 
functional form of the covariance between the number of eigenvalues in overlapping 
intervals of equal size. Following up on 07], where Odlyko studies the difference 
between the spacing distribution for the Riemann zeros and the corresponding lim¬ 
iting random matrix distribution and writes: “Clearly there is structure in this 
difference graph, and the challenge is to understand where it comes from”, Bo- 
gomolny et al. [8] considered finite size corrections to the spacing distribution for 
random unitary matrices. While previous studies in random matrix theory gave 
analytic results relating to the value distribution of the characteristic polynomial 
[T], there is no existing literature on the analytic calculations of the leading finite 
size correction to the spacing distribution for random unitary matrices. Here we are 
referring to the point-wise value of the expected probability density function; there 
is an existing literature relating to the decay as a function of N of the Kolmogorov 
distance between the empirical cumulative spacing distribution and its expected 
large N form (SSJ |32] . In |S] the corrections were computed using an extrapola¬ 
tion procedure of a determinant expression valid for finite N. A primary objective 
of the present paper is to provide an analytic characterisation of the leading finite 
size correction to the spacing distribution. 

Variants of the nearest neighbour spacing distribution also provide relevant sta¬ 
tistical quantities. For example, one could consider spacing distributions between 
zeros k apart (see e.g. m §8.1]), or the distribution of the closest of the left and 
right neighbours [23]. Another variant is to consider spacing distributions resulting 
from first thinning the data set by the process of deleting each member indepen¬ 
dently with probability (1 — ^), 0 < ^ < 1, as first considered in random matrix 
theory by Bohigas and Pato m- In both cases the finite size correction to the 
large N form of the corresponding quantity for the eigenvalues of random unitary 
matrices can be computed analytically. This is of concrete consequence as we are 
fortunate enough to have had A. Odlyzko provide us with a data set extending that 
reported in m- Specifically, this data set begins with zero number 

100,000,000,000,000,985,531,550 « 102^ 
and this occurs at the point s = 1/2 + iE in the complex s-plane with E equal to 
13066434408793621120027.3961465854 « 1.30664344 x 10^^ 


The data set contains a list of slightly more than 10® of each of the subsequent 
Riemann zeros. Thus we are able to compare the analytic forms against the finite 
corrections exhibited by the Riemann zeros, uninhibited by sampling error. 

We begin in Section]^ by considering the two-point correlation. Bogomolny and 
Keating m have derived an explicit formula for this quantity in the case of the 
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Riemann zeros, in the regime that E is large but hnite. Subsequently Bogomolny 
et al. [5] showed how the result of m could be expanded for large E to obtain 
the correction term to the random matrix result (1.1). Moreover, upon (partial) 
resummation, the correction term was found to be identical in its functional form 
to the correction term of the scaled large N expansion of the two-point correlation 
for the eigenvalues of random unitary matrices chosen with Haar measure. After 
reviewing this result, we compare the empirical form of the two-point correlation as 
computed from Odlykzko’s data set, and for the data set thinned with parameter 
f = 0.6, against the theoretical prediction. Very good agreement is found for 
distances up to approximately one and a half times the average spacing between 
zeros, after which systematic deviation is observed. 

Higher order correlations for random unitary matrices are fully determined by 
the same function — see (2.1) below — as determines the two-point correlation. 
Assuming the same is true for the correlations of the Riemann zeros, it follows that 
the coincidence between the correction term to the leading large E form of the two- 
point correlation for the Riemann zeros and the correction term to the leading large 
N form of the two-point correlation for the eigenvalues of random unitary matrices 
persists to all higher order correlations. While higher order correlations cannot 
be measured empirically from the Riemann zeros data, this hypothesis becomes 
predictive as the coincidence must carry over to any distribution function which can 
be written in terms of the correlation functions, for example the nearest neighbour 
spacing distribution [8] . In Section 3 we make two main contributions to the study 
of this theme. One is to test the hypothesis upon thinning of the data set. The 
other is to provide an analytic determination of the correction term in the random 
matrix case, using certain Painleve transcendents. 

The correction term to the nearest neighbour spacing is oscillatory, and becomes 
more so as the thinning parameter ^ is decreased from 1. In Sections and 
we discuss this effect in the context of the large distance form of the Painleve 
transcendents. 


2. Two-point correlation 

Correlation functions are fundamental to the theoretical description of general 
point processes. For definiteness, and in keeping with the setting of the Riemann 
zeros, we will specify that the point process is defined on a line. Starting with 
/9(i)(a;) as the density at point x, the fc-point correlations ...,ccfc) can be 

defined inductively by the requirement that the ratio 

P{k){p^l^ ^2 ; • ■ • 7 ^/c) 

P(fc—1) (^1; ^2; • ■ ■ j l) 

corresponds to the density at the point Xk, given that there are particles (or zeros, 
or eigenvalues etc.) at locations Xi,... ,Xk-i- Suppose the particle density /5(i)(x) 
is identically constant so that the point process is translationally invariant, and 
furthermore take the constant to be unity. We then have that the two-point corre¬ 
lation /3(2 )(xi,X 2) is just the density at X 2 given that there is a particle at xi. As 
such, in this circumstance, p(i)(xi,X 2 ) can be empirically determined from a sin¬ 
gle data set. To see this note that due to translational invariance, p( 2 )(xi,xi -I- s) 
depends only on s, allowing xi to be averaged over to generate the necessary data 
for a statistical determination. 
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For the eigenvalues of iV x unitary matrices chosen with Haar 

measure, the fc-point correlation is given by the k x k determinant (see e.g. |271 
§5.5.2]) 

P(k)ixi,. ..,Xk) = det[KNix^,Xj)]^j^i^„,^k 

where Xj = N9j/2'K (j = 1,... ,7V) and Kjs[{x,y) — referred to as the correlation 
kernel — is given by 

sin7r(a; — y) 


KNix,y) = 


N sm{TT{x — y)/N) 


( 2 . 1 ) 


smTTS 


In particular the two-point correlation is given by 

P( 2 )ix, a;-|-s) = l— 1 

^ \NsimTs/N^ 

which one sees tends to ([H in the limit TV —>■ oo. We are particularly interested 
in the leading correction term to this asymptotic form. As reported in [5], an 
elementary calculation gives 

, , , /sin7rs\2 1 „ 1 \ 

P( 2 )(x,a: + s) = l-(^^j sin Ts + 0(^J. (2.2) 

The situation with the Riemann zeros is more complicated. For a start, since 
the Riemann zeros are a deterministic sequence, a statistical characterisation only 
results after an averaging over a suitable interval of zeros about the point x of 
interest along the critical line, and with the use of a test function. And, as men¬ 
tioned in the Introduction, it is only for a restricted class of test functions — those 
whose Fourier transform have support on jfcj < 2?! — that rigorous analysis of the 
correlations has been possible. Even then the theorems obtained are statements 
about the leading asymptotic form only. 

The discussion in the Introduction mentioned that a non-rigorous, but predictive 
analysis based on an analogy with a chaotic quantum system, together with the 
use of a conjecture of Hardy and Littlewood for the pair correlation function of 


the primes, allows for further progress m- One consequence is that ( |I.2[ ) can be 
derived as the two-point correlation function without the assumption of a restricted 
class of test functions. Again, this is a statement about the limiting asymptotic 
form. For present purposes, the most relevant feature of the results of m is that 
they allow for the determination of the leading correction term to the limiting 
asymptotic form. The necessary working is given in [S], where it was shown that 
at position ^ -|- iE, with E ^ 1, along the critical line, with the local density 


normalised to unity, and with p given by the leading term in (I.I) the smoothed 
two-point correlation has the large E expansion 


R ^ , /sinTTSN" 

P( 2 )iE,E + s) = 1- 


A 

■k'^P^ 


sm TTs — 


Q 


rS sin 27rs 


0(i). (2.3) 


2'k'^P^ 

The constants A and Q can be expressed in terms of convergent expressions in¬ 
volving primes, which when evaluated give the numerical values A = I.573I4... and 
(5/A = 1.4720.... Moreover, with 


a = I 


C 


log(A/27r) ’ 


C = 


Q 

A’ 
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it was observed that (2.3) can be rewritten 

, 2 


R/ N , /sinTTSN^ A .9, , ^/i\ 

+ = -^sin (7ras) + 0(^j. 




1 


(2.4) 


Comparison with ( |2.2[ ) shows that the correction terms agree subject to relating 
the matrix size N and the distance along the critical line E according to [8] 


yiM V27r/ 

provided too that the scaled distance s is further rescaled 


(2.5) 


s i-A as. (2.6) 

In the theory of point processes (see e.g. [36]) a thinning operation, whereby each 
member is deleted independently with probability (1 — ^), is used to create a family 
of point processes from a single parent process. The effect on the corresponding 
fc-point correlation P(k){xi,.. .,Xk',^) is very simple — the ^ dependence scales to 
give 

Pik){xi,...,xk;0 = ^’'P{k){xi,...,Xk;l). (2.7) 

Thus introducing the variables yj = ^Xj {j = 1,2,...) one sees that 

P{k) {xi,..., Xk\Ofi'Xi ■■■dxk = P(k){yili, ■■■,yk/^; l)dyi • • ■ dyk- (2.8) 

The variables {yj} correspond to a rescaling so that the density of the original point 
process remains unchanged, assuming translational invariance. 


2.1. Numerical results. To determine the two-point correlation P( 2 )(?/i/^; 2 / 2 /^; 1) 
empirically from Odlyzko’s data set we must first rescale the data by so that 
the rescaled local mean density is unity; we then delete each zero with probability 
1 — Next we sample this sequence by empirically computing the zero density 
of the sub-system formed by each zero in turn as the left boundary, and a hxed 
number (say 50) of its neighbours. Figure displays the results of this empirical 
determination of (2.8) in both the variables {xj} (on the left) and {yj} (on the 


right). To the eye, the resulting graphical forms are identical to the leading order 
random matrix prediction (1.2). 


Following [^ our main point of interest is in the functional form that results 
after subtracting the leading random matrix prediction from the empirical two- 
point function for Odlyzko’s data set. 


.g,(E.^+»)-(i-(5;^) 


(2.9) 


The theory of [5|, as reproduced in (2.4) above, predicts that this quantity is, to 


leading order in 1/p, itself of random matrix origin. As a test of its validity. Figure 
id isplays a comparison of the empirical data and the theoretical predictions. Very 
good agreement is found for distances up to approximately one and a half times the 
average spacing between zeros, and furthermore the period of the oscillations are 
well matched for all distances in the display. On the other hand, as the distances 
increase, there is a systematic discrepancy at a quantitative level. Anticipation 


of the expanded form (2.4) not being in agreement for increasing distances comes 


when one graphically compares the corrections in (2.3) and (2.4) — recalling, in 


particular, that theoretically (2.3) and (2.4) agree to leading order. Thus one finds 
that their mismatch increases with increasing distance. On the other hand, the 
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functional form from m used to obtain ( |2.3| ) can itself be numerically calculated. 
Its graphical form has for ^ = 1 been compared against the two-point correlation of 
the same sequence of Riemann zeros as used here [SJ Figures 12-17] or [71 Figures 
3-7] and quantitative agreement is found for all distances displayed. 




Figure 1 . On the left is a plot of the raw 2-pt correlation function 
in the Riemann zero data, where each point has been deleted with 
probability (1 — ^). The figure on the right contains the same 
data, which has been rescaled according to (2.81, showing the self¬ 
similarity property for varying 



Figure 2. Comparison of the 1/p^ term in ( 2.4 ) [solid curve] (with 
s —> s/^) against (2.9) [crosses] (again with s —>■ s/^). 


3. Spacing distributions 

3.1. Fredholm determinant formulae. Although the fc-point correlations can¬ 
not be empirically computed from the data beyond the case fc = 2, certain function¬ 
als of the correlations can be so computed. A specific example is the n-th nearest 
neighbour spacing distribution p{n; s), n = 0,1, 2,.... In terms of the correlations, 
and assuming a translationally invariant system with constant density p, we have 
in the case n = 0, 

p(0; s) = Q) 1 ^ f dxi--- [ dxj p(^j+2){s,0,xi,...,Xj), (3.1) 

P ^ j=i JO Jo 

and similar formulae for general n; see e.g. [271 §8.1 & §9.1]. 

The calculation of the functional form of the scaled limit of p(0; s) for A^ x A^ 
unitary matrices U(N) chosen with Haar measure, or equivalently bulk scaled CUE 
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matrices, is a celebrated problem in random matrix theory. First note that a 
Fredholm determinant can be expanded as a series involving multiple integrals 
with integrands given as determinants of the integral kernel. With this expansion 
one can show from the form of the scaled fc-point correlation functions (1.3) that 

(3.2) 

where Kg is the integral operator supported on (0, s) with kernel 


p(0;^ fiet(I -: 


K(x,y) = 


SinTT 


{x - y) 


tt{x - y) 

noting that the constant density p is unity. The required working can be found, 
for example, in §9-1]. This Fredholm determinant formula is equivalent to the 
expression 

j 2 °° 

1=0 

where 1 > Ai(s) > A 2 (s) > • • • > 0 are the eigenvalues of Kg. Such a result, for 
the bulk scaled limit of GOE matrices, was first obtained by Gaudin [53], who 
also provided a computational scheme by relating the eigenvalues {Ai(s)} to the 
eigenvalues of the second order differential operator having the prolate spheroidal 
functions as eigenfunctions (see e.g. [23 §9.6.1]). Only recently has it been shown 
that there are numerical schemes based directly on (3.2) that exhibit exponentially 
fast convergence to the limiting value mull. 

Nearly two decades after the work of Gaudin, the Kyoto school of Jimbo et 
al. [33 expressed (3.2) as the so-called t- function for a particular Painleve V system. 
Explicitly, it was shown that 

a(0)(t;e) 


det(I — ^Kg) = exp 

where a satisfies the differential equation 


dt, 


4:{ta' — a){ta' — cr -I- {a')^) = 0 


with small t boundary conditions 

^( 0 ) 


TT TT^ 




(3.3) 


(3.4) 


(3.5) 


Note that the parameter ^ introduced in (3.3) only enters in the characterisation 
through the boundary condition. Eor background theory relating to the Painleve 
equations as they occur in random matrix theory we refer to [23 Gh. 8]. 

Although (3.2) only requires the case ^ = 1, the Eredholm expansion formula 
used to derive (3.2) generalises to read 

i-0^ 


J"_ ^'\J rS nS 

det(I - ^Kg) = ^ —— / dxi--- dxj pQ){xi,... ,Xj), 

. n 3- Jo Jo 


j=o 


(3.6) 


(cf. eq. (3.1)). In the case ^ = 1, this formula, via the inclusion/exclusion principle, 
shows det(l — Kg) can be interpreted as the probability, E(0; s), say, that there are 
no zeros in an interval of length s of the original data set. Thus we have 

E{0; s) = det(I — Kg). 


(3.7) 
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In the case ofO<f < 1, (2.7) together with the inclusion/exclusion principle shows 
(3.6) is equal to the probability that there are no zeros in an interval of length s in 
the data set corresponding to a thinning of the original data set by deleting each 
zero with probability (1 — ^) — this is the procedure alluded to in the Introduction. 
With p(0; s; denoting the corresponding distribution of the nearest neighbour 
spacing, (|3.2[) generalises to 


1 (f 

p{0; s; 0 = ^ ^ det (I - , 


(3.8) 


where the the constant density p = ^ in this decimated system. We remark that the 
notation p((0, s);^) has been used in [571 §8-1] to denote the generating function 
for the conditioned spacing distributions {p(n; s)}„=o,i,...j corresponding to the 
conditioned gap probabilities in the next paragraph. These quantities are related 

by 


1 


p(0;s;^) =p((0,s);O. 


(3.9) 


As just mentioned, another interpretation shows itself via the formula (see e.g. da 
eq. (8.1) together with (9.15)]) 


E{m] s) = 


(_l)m gm 

ml 


■ det (I - ^Ks) 




(3.10) 


where E(rn]s) denotes the probability of an interval s in the original data set 
containing exactly m zeros. This is referred to as a (conditioned) gap probability. 
Note that (3.10) reduces to (3.7) in the case m = 0. 

Both the expressions (3.2) and (3.3) have analogues for the eigenvalues oiU{N) 
itself, rather than their scaled limit. Let Kf denote the integral operator supported 
on (0,s) with kernel Kjs[{x,y) as specified by (2.1). The analogue of (3.8) is then 
the structurally identical formula 


1 d- 


p^(0;s;C) = ^^det(I-eKf). 


(3.11) 


Here p^(0;s;^) denotes the nearest neighbour spacing distribution for matrices 
from U(N) chosen with Haar measure, with s the difference in neighbouring angles. 
The T-function formula (3.3) generalises to |32l eq. (1.33)] 


I-ks/N 

det(I — ^Kf^) = exp — y U{cot (j>-, d(j>j, 
where U satishes the (fPVI equation (see e.g. m eq- (8.21)]) 

4 

m'(( 1 + s^)u")^ + 4(zi'(u - su') + iViV 2 V 3 V 4 )'^ + 4 (u' + vl) = 0, 


(3.12) 


(3.13) 


k=l 


with u{s) = U{s;(,) and parameters 

Vl = V2 = V3 = 0 , 


V4 = N, 


subject to the boundary condition u{s) ~ £,N/tt as s —>■ oo. 
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3.2. Finite N correction for the nearest neighbour spacing. Our interest is 
in the leading correction term to the large N form of (3.11). Since, from (13, 


simr{x-y) tt^x - y) 
KN(x,y) = -^^-h 


/ N I ■ sinTrta; — n) + 0( —7 I 

n{x-y) 6iV2 ^ \N'^J 

we see that this correction term is of order 1/iV^. Specifically 


det(I-CIKf) = det(I-CIK,) 


( 3 . 14 ) 

where Lg is the integral operator on (0, s) with kernel (7r(x — j/)/6) sin7r(a: — y). 
We remark that with R{x,y) denoting the resolvent kernel — that is the kernel 
supported on (0, s) of the integral operator ^IKs(I — — straightforward 

manipulation shows 

Tr(^(I-^Kg)"^Lg^ ^( 2 /-a;)- 


The expansion (3.14) tells us that in relation to the representation (3.12) we 


should change variables and write 


+ 0 




V 


(3.15) 


This implies the expansion 
det(I — ) = exp 


r(0)(X) 


dX 






dX + O —J 


(3.16) 


Both and also depend on ^ (which, as mentioned below (3.5), enters 
through the boundary conditions), but for notational convenience this has been 
suppressed. As noted above, satishes the particular Painleve V equation in 
sigma form (3.4) with boundary condition (3.5). By changing variables s = cot X/N 
in ( |3.13[ ) and introducing the expansion ( |3.15[ ), the equation (3.41 can be reproduced 
(a fact already known from |31jl. and moreover a linear differential equation for 
(T*'^^(A) with coefficients involving can be obtained. 


Proposition 3.1. With a^^'^{X), a^^\X) related to U{s;^) in (3.12) by (3.15) 
we h ave that ct^°^(A) satisfies the particular Painleve V equation in sigma form 
(3.4) with boundary condition \3.t^ , while ct*^^^(A) satisfies the second order, linear 
differential equation 


A{s)y''{s) + B{.s)y\s) + C{s)y{s) = Dfs), (3.17) 

where, with cr{s) = (T*-°^(s), 

A{s) = 2 s^(t"(s), 

B{s) = —8a'{s)a{s) + 12s(ct'(s))^ + 8s^s(t'(s) — cr(s)^ , 

C{s) = -A{a'{s)f - 8(^sa'{s) - cr(s)), 

D{s) = -^sV"(s) ^cr(s) - sa'{s) - 

- ^(scr'(s) - (T(s))(^3(cr(s))^ + 2scr(s)(s - cr'(s)) - 2 s^(t'(s)(s + ^'(s))^ 
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The corresponding s —>■ 0“'' boundary condition is 

, „5 5^' 




(3.18) 


Proof. We see there are three main terms in (3.13). Changing variables s = 
cotX/N and introducing the expansion (|3.15[) in the first gives 


(^{l + s^)U"is;0y =N^ 




dr (7 




J 


+ AX^ 


/I Xda^^^X) a(°\X) X^ {X) 


dX2 


dX2 


dX 


6 dX^ 


up to terms of order 1/N'^. Similarly expanding the other two main terms and 
then equating terms of order N'^ gives (3.4), while equating terms independent of 
N gives (3.17). 


In relation to the boundary conditions, for U{s;f ) in (3.12) we know from [23 
eq. (8.78)], further extended to the next order, that for large s 


U (s; —c H-1—y + ' 


, 9c^ -N^ -2 
9s^ 


. 36c2 - 5N^ - 19 

36?^ 


2 46 - 375c2+450c4 + iV2(40 - 75c2) + 4Af4 

+ ■= -455?^-^ 


where c = ^. Substituting in (3.15) and performing appropriate additional ex¬ 
pansions we read off (3.18). □ 

While we know of no other analytic formulae for the correction term to the spac¬ 
ing distribution for random matrix ensembles in the bulk, at an edge — specifically 
the soft edge — this task seems to have been first taken up by |2D]. Aspects of this 
same problem at the hard edge are discussed in the recent works (23 [HI EH] . 

Let {x)[N~P] denote the coefficient of N~p in the large N expansion of {x). 
Then we have from (3.16) that 


det(I-eiKf)[fV-2] = (3.19) 

One use of Proposition is to be able to deduce the explicit form of the small and 
large distance expansions of (|3.19 ). We will consider first the small s expansion. 


Corollary 3.2. We have 


f2„2 


71 ® 677 '* 225772 J 

121^* 


24 V3t2 


24f 

'7r4 


s'‘ -f 


(5e e 

y36772 77' 


7^ _ 7C^ \ 
77^ ' 3677® 67577® J 


' 778 9776 

+ 0(si°). 


-f 




810077* 2205772 


- 




73C 


761^® \ 


2 477'^ 360077® 158760077® / 


(3.20) 
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+ , 7(-15C® + 2CV^)s^ 

' “ n_2 


SGtt® 


QOtt^ ^ 5407 r® 

(1260C® - 2037r2C'‘ + 127r^C^) s® (9450^'^ - ITSStt^^® + SStt^^®) s® 
56707r6 " 

+ 0(^1®), 


378007r7 


(3.21) 


det(I-C]Kf)[l] = 1-^5 + 


36 




c2„6 
? ^ 


? tt 


s” + " • s® + O (s®°) , 

675 17640 291600 ^ ’ 

(3.22) 


t2_2„4 c2_4„6 c 2 ^ 6„8 <r3_6 „9 

det(I - ^Kf )[Ar-®] = ^ + O (s“) , (3.23) 

J 36 270 3780 48600 '' '' ^ ^ 


27r 

iV 


p^\0-,27rs/N;^)[l] = 




3 


0-7^ 

-p^(0;2WYO[iV-^]=- 3 


45 

Ctt^s^ 

9 


315 

2^^656 

135 


4050 


o 


(3.24) 


675 


+ 0 (/). 


(3.25) 


Proof. The differential equation (3.4) admits a unique power series expansion 


about the origin with the first two terms given by (3.5). Expanding to higher order 


and solving for the unspecified coefficients gives (3.20). Substituting the expansion 


in (3.3) gives (3.22) for the small s expansion of the Fredholm determinant; this 


series is also reported in m eq. (8.114)]. If we substitute ( |3.20[ ) instead in the 
differential equation (3.17) and make use of the boundary condition (3.18) we find 


a unique power series expansion solution is generated which upon solving for the 


unspecified coefficients gives (3.21). Substituting (3.20) and (3.21) in the order 


uniform density P(i)(s) = we see (3.24) and (3.25 
(Iff^ and plt|. 


l/iV2 correction term ( |3.19 ) gives (3.23). Recalling (3.11) and that we have a 


follow immediately from 

□ 


A consistency check can be placed on the abo ve ex pansions. For this we note 
that it is also possible to use the characterisation (3.12) of det(I — to deduce 
the A^-dependent small s expansion [33 eq. (8.79) 

(l-l/fV2)^2^2^4 


det(I-CKf) = 1-Cs + 


(l-l/Af2)(2-3/fV2) 


36 


1350 




(l_l/iV2)(l-2/iV2)(3-5/iV2) 


52920 


^2^t>s» + 0(s®). (3.26) 


Expanding the RHS of (3.26) for large N we obtain agreement with (3.22) for the 
leading form, and agreement with (3.23) for the next leading order 1/iV^term. 


A noteworthy feature of (3.24) and (3.25) is the weak dependence on the dilution 


parameter which only changes from ^ to at order s^. One understanding of this 


result relates to the interpretation of (3.8), upon division by as the generating 
function for {p (ri,; s) }„-o 1 ; r ecall (3.9). Then, in analogy with (3.10), p(l;s) is 

obtained from (3.24) and (3.25) by applying —^ and setting ^ = 1. It follows that 
p(l;s) has leading small s behaviour proportional to s^, which is a known result 
eq. (8.115)]. 

We turn our attention now to the behaviours for large s. The asymptotics 
of spacing distributions for random matrix ensembles in this regime have been 
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reviewed in the recent work [28]. The case ^ = 1 must be distinguished from 

0 < e < 1. 


Corollary 3.3. Suppose 0 < e < 1. s 
cro(s/2) = cr(°)(s) in our notation] 


a^°\s) = —ks + + O , k 


have [231 eq. (1.16), with 
^ (3.27) 


=-log(l - C) 

TT 


and thus 

Also 


det(I-CKf)[l] ^ 7l(0s--e 


j2 —kirs 


cr^^^(s) = — 

and this together with 28) implies 


k^s^ 


0(s) 


det{l-^K")[N-^]= Ai^s'^ l^e 


k^ j2^—k'KS I _ ^ 

12 


(tts)^ + 0(s) 


Suppose instead ^ = 1. In this case, for s —>■ oo we have 


r(0) 


1 

(s) = -^-7 


0(1 


which is equivalent to the expansion 

det(I-Kf)[l] 

We also have 


A{1) 


„-{TTsf/8 


.(1) 


/ \ 


s 


which when combined with (3.32) is equivalent to 


det(I-CKf)[lV-^]=ll(l) 
Thus, for large s 


(tts)!/^ 


(7rs)i/4 

0 ( 1 ), 

(tts)"^ ("^s)^ 


192 


96 


0(5) 


27r 


p"(0;2Ts/iV;O[l] 


A(l) 


(k'Kfs^''l'^e-^^\ 0 < ^ < 1 , 


IL^c7/4„-(^s)V8 


16 


?=1, 


^TT 

-p^(0;27rs/lV;^)[Ar- _ 
iV s—>-oo 


r-2i 


0<C<1, 

-ll(l)l^s23/4e-(-)V8, 


e = l. 


(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3.35) 

(3.36) 


Proof. As noted, the expansion (3.27) can be read off from [43], and this substi¬ 
tuted in (3.3) gives (3.28). The value of A(^) is given in [T3 eq. (1.14)], although 
its derivation requires different methods [131 m . Substituting ( |3.27| ) in (|3.17| ) and 
solving for large s gives (3.29), and this together with knowledge of (3.28) gives 
( |3(30| ). 

From the definition of k, we see that the behaviour (3.27) breaks down when 
^ = 1. In that case we have, instead of (3.27), the large s expansion (3.31). This 
follows from ( |3.3| ) and knowledge of the leading two terms of the large s form of 
det(I —Kg) |24L 154] as given by (3.32). Substituting (3.31) in (3.17) and solving for 






























14 


PETER J. FORRESTER AND ANTHONY MAYS 


large s gives (3.33). Then substituting this in (3.19) and using (3.32), we see that 
the large s form of the 1/N'^ correction term (3.19) for ^ = 1 is given by (3.34). 


To obtain the behaviour of the spacing distribution we use (3.28) 


3.28 

I, 

3.30 

I, ( 

3.32 


and (3.34) in conjunction with (3.11) and so deduce (3.35) and (3.36). 


□ 


Numerical methods to be discussed in Section [fallow the large s forms of (t*'°^(s) 
and (s) to be compared against their computed values from the differential 
equations; see Figures [SJQ A prominent feature seen upon differentiating the 
transcendents is that for 0 < ^ < 1 higher order terms in the large s expansion are 
oscillatory, as known from the analytic work of McCoy and Tang [33] in the case 
of cr*^°^(s). Since there is no such effect for ^ = 1, we display only the transcendent 
and its asymptotic form, and not the derivative. The asymptotic form in this case 
is functionally distinct from that for 0 < ^ < 1 — note that the quantity k in the 
latter actually diverges at ^ = 1. 




Figure 3. ^ = 0.6. The black curves are the solutions of the DEs 


( |3.4| ) [left] and (3.17) [right] using a sequence of Taylor series of de¬ 
gree 12. The dashed red curves are the asymptotic approximations 
(3.27) [left] and (3.29) [right]. The inset is the same plot zoomed-in 
on the origin. 



corresponding curves in Figure]^ 

A variation on the process of deleting each zero with the (constant) probability 
1 — ^ is to delete with probability 1 — ^(cc), with x = 0 corresponding to some pre¬ 
determined origin in the data. Since upon the scaling (O) the original data set is 
translationally invariant, an ensemble can be generated by varying the origin. The 
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Figure 5. ^ = 1. The black curves are the solutions of the DEs 
( |3.4| ) [left] and ( |3.17 1 [right] using a sequence of Taylor series of 
degree 35. The dashed red curves are the asymptotic approxima¬ 
tions (3.31) [left] and (3.33) [right], 
zoomed-in on the origin. 


The inset is the same plot 


corresponding A:-point correlation, (cci,..., Xa,; ^(a;)) say, has the factorisation 
property 


P(k){xi,... ,Xk]i{x)) = p(k){xi,...,Xk\l) 


(cf. (2.7)). Thus in this circumstance the correlation kernel in (1.3) should be 
replaced by 


{i{xj)£,{xi)) 


i/2Sin7r(a:j — xi) 
Tr{xj - xi) 


(3.37) 


Denoting the integral operator supported on (0,s) with kernel ( |3.37 1 by IKs[^], 
the probability that an interval of size s from the origin has no zeros is equal to 
det(l — Ks[^]). Subject to ^(x) having the functional form ^(x/s) = u(y), y = x/s 
with v{y) analytic on [0,1], then (for large s) the expansion (|3.28| generalises to 
[m Theorem 2.1] 

det(l - K4^]) - 


where 

1 /"i 1 

ky = -- \og{l - v{x))dx, ay = - \og{l-v{y)), 

71 7o tt 

and for a certain explicit C'[^]. Unfortunately the requirement that ^ scale with s 
prohibits probing this process in the Riemann zeros data. This is similarly true of 
the choice = 1 — recently studied for large s in [18] . 


3.3. Alternative characterisation of finite N correction for nearest neigh¬ 
bour spacing. We see from (3.11) that to obtain the order 1/A^ correction to 
the limiting spacing distribution from knowledge of the 1 /7V^ correction to the gap 
probability, the second derivative with respect to s of the latter must be taken. In 
fact, as first shown in |30| and further refined in m, the second derivative can 
be incorporated within the Painleve theory. Specifically, we have from m (again 
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recalling (3.9)) that 

P{0;s;^) =—^exp J 

where satisfies the particular, modified crPV equation 

^su"(s)^ + ^su'(s) — u(s)^ ^su'(s) — u(s) — 4 + 4 (m'(s))^^ — 16^u'(s)^ =0 

(3.38) 

subject to the small s boundary condition 

;(s® + 0(/)^ 


Moreover, the result m eq. (5.16)] also tells us that for finite N 


(3.39) 


^p^(0;27rs/iV;^) = - 1) sin^ ^ exp |- 


s/N 


y(cot(^;^) dij) 


where V satisfies the ctPVI equation (3.13) with u{s) = F(s;^) and parameters 

vi = V4 = 0, V2 = —N, Vs = —2. (3.40) 


Analogous to (3.15), in keeping with the correction to the large N form being of 

1 u(i)(X) 


order l/N'^, we expand 


-l^VicotX/{2Ny,0 = -^ 


Af2 X 




so that 




O-rr 

—p^(0; 2tts/N] exp 


X 1- 


^27rs 


TT^S^ 


j(o)(X) 

1 f 


dX 


1 


N2 3jV2 iV2 


X 


dX + O 


1 

iV4 


(3.41) 


We know that m^^X) satisfies the nonlinear e quati on (3.38) subject to the bound¬ 
ary condition (3.39). Analogous to Proposition 3.1 it is possible to specify u^^^(X) 
as a second order linear differential equation with coefficients involving u^^\X). 

Proposition 3.4. IPe have that u^^^(X) satisfies the second order, linear differen¬ 
tial equation 

A{s)y"{s) + B{s)y'{s) + C{s)y{s) = D{s), (3.42) 

where, with u{s) = u^^\s), 

A(s) = 8s^it"(s), 

B{s) = 8 (6s(u'(s))^ -I- s^m'(s) — 2s — su{s) — 16u'(s) — 4u(s)u'(s)) , 

C{s) = 8 (2 -I- u{s) — su'{s) — 2{u (s))"^) , 

D{s) = ^ s^u"(s) ^s^u"(s) -|- 2su'{s) — 2u{s)^ s'^u'{s)^ 4s^u'(s)^ 

— 2s^u{s)u'(s) — 2s^u'{s) 16s^u\s)^ s^u{s)^ 2s^u{s) 

— 10su(s)^u'(s) — 64su(s)it'(s) — 96su'(s) -h 6u(s)^ -h 48u(s)^ -h 96u(s) 


















FINITE SIZE CORRECTIONS AND ODLYZKO’S DATA SET FOR THE RIEMANN ZEROS 17 


The equation must be solved subject to the s 


,(i) 




13 


O'*" boundary condition 


0(s‘>). 


(3.43) 


15 6300 172877 

Proof. We apply the same technique as in Proposition [3^ to the trPVl equation 


(3.13) with parameters (3.40); that is we change variables s = cot X/{2N) and 
we find the terms of order give ( |3.38 ), while the terms of order 1 give 
For the boundary condition on we have (3.39), which we extend to degree 


6 using (3.38). Then we combine this series with (3.42) to obtain the boundary 
condition on (s). □ 


One check on the consistency of Proposition 3.4 is to verify that substitution 


of (3.43) and the first few terms of the power series solution of (3.38) into (3.41), 


reproduces the expansion (3.25). We find that this is indeed the case. 

The required consistency between (3.41), and (3.16) combined with ( |3.11[ ), to¬ 
gether with knowledge of the large s asymptotic forms in Corollary |3.3| give infor¬ 
mation on the large s behaviour of and In relation to for 

0 < 5 < 1, we have [331 eq. (1.26), with d 2 {s/f) = u^^^(s) in our notation] 

= + + (3.44) 


(up to a correction bo(n) —bo(n) in |43l eq. (1.13)j) with k as in (3.27). From 
[43] we know that the 0(l/s) term in ( |3.44 ) is oscillatory. This is clearly displayed 
in Figure [^ where we compare a numerical solution of the differential equation 
(3.38), with [331 eq. (1.26)] as given in (3.44) with the oscillatory term 0(l/s) 
included. These oscillatory terms have concrete significance for the deduct ion of 


the functional form of u^^\s, in the large s limit, as if we are to substitute (3.44) 


for this purpose in the differential equation (3.42), we find a result inconsistent with 


(3.36). Consistency with the latter and (3.29) requires 


,(i) 


(5,0 = y + 0(s). 


(3.45) 


Note that this leading term is independent of From Figure [^ we also see that 
there are oscillatory sub-leading terms in (3.45); the graphs of the derivatives in 
Figure [^indicate that these oscillations begin with the 0(s) term. 

For ^ = 1, we have consistency with (3.35) when 


' ^ 16 4 


O ( - 
s 


(3.46) 


and substitution of this in (3.42) gives 


,(i) 


(s,l) = -I 


0 ( 1 ); 


(3.47) 


43s^ 

768 J92 

see Figure for comparisons of these asymptotic forms against numerically gener¬ 
ated solutions of the corresponding differential equations. 

4. Numerical power series solution of the differential equations 
AND evaluation OF THE SPACING DISTRIBUTIONS 


The utility of Propositions |3.1| and |3.4| is that we can use a power series method 
together with computer packages to numerically compute to high accuracy the next- 


to-leading order term in the nearest neighbour spacing distribution (3.8); this can be 
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Figure 6. ^ = 0.6. The black curves are numerical solutions of 
the DEs (3.38) [left] and (3.42) [right] using a sequence of Taylor 
series of degree 12. The dashed red curves are the asymptotic 
approximations (3.44) [left] and ( 3.45| ) [right]. The dashed blue 
curve in the graph on the left is a plot of the three leading terms 
in [331 eq. (1-26)] with n = 2,9 = 0,k = ^ log(l — therein. 




Figure 7. ^ — 0.6. These are the derivatives w.r.t. s of the 
corresponding curves in Figure [^ 



Figure 8. ^ = 1. The black curves are numerical solutions of 
the DEs (3.38) [left] and (3.42) [right] using a sequence of Taylor 
series of degree 35. The dashed red curves are the asymptotic 
approximations (3.46) [left] and ( |3.47| ) [right]. The insets are the 
same plots zoomed-in on the origin. 


done directly by computing the terms in (3.41) or via the second derivative of (3.19). 
Such power series methods were introduced as a technique to compute Painleve 
transcendents associated with spacing distributions in random matrix theory in 
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[49]. This was in relation to the soft edge. Subsequently the same numerical 
method was used for the computation of bulk spacings [13 §8.3.4]. The first point 
is to note that in generating a power series solution of the crPV equation (3.4) 
about the origin, one finds a radius of convergence of approximately 8.5 (for both 
^ = 1 and ^ = 0.6). Inside the radius of convergence, this power series can be used 
to accurately evaluate the transcendent and its derivative, and from this data a 
new power series can be computed and the procedure iterated to cover the interval 
from s = 0 to beyond s = 20. With (T*-°^(s) so determined, iterative power series 
solutions of the differential equation in Proposition 3.1 can be obtained. It is 
through this procedure, and its analogue starting with the differential equation 
(3.38) and proceeding to the differential equation of Proposition 3.4 that the graphs 
of (s), (s), (s) displayed in Figures [5||3 have been generated. 

We remark that use of a power series method to generate solutions over a large 
interval seems necessary in the cases 0 < ^ < 1. In particular, if in these cases 
one tries to use a computer algebra package to solve the differential equation (3.4) 
subject to initial conditions for a^^\xo) and cr^^^'(xo) with xq small, as computed 
from the boundary condition ( |3.5| ) (with the latter further extended in accuracy 
as given in (3.20l), it is found that as x is increased from xq, the DE solver gives 
incorrect values, and furthermore soon diverges to a spurious pole. _ 

Substitution of the piecewise functions for tT(°^(s) and cr(^)(s) into (3.16) gives us 
an approximation of det(I — up to the first order correction, and substitution 

into ( 3.19| ) isolates the correction term. Taking the second derivative with respect to 
s we obtain Figures!^ and [Toj respectively. Note that to obtain the comparison with 
the Riemann zero data in these figures, we have used the scalings (2.5) and (2.6) 
from [Sj. Also in Figures [^ and [To] we plot an extrapolation of a sequence of finite 
calculations ( 3.11| ), similar to that in [5]: by calculating (3.11), with kernel (2.1), 
for 20 values of N between 100 and 138 we extrapolate the limiting value of the 
nearest neighbour spacing and the next-to-leading order correction at each point. 
To graphical accuracy the extrapolated values are identical to those computed from 
the differential equation. Note that when ^ = 0.6 the correction term shows more 
pronounced oscillations for finite s values than its ^ = 1 counterpart, a feature that 
seems to be be driven by oscillations in the functional forms for the corresponding 
Painleve transcendents. 

Most strikingly, the predicted functional forms from random matrix theory show 
accurate agreement with the Riemann zero data both for ^ = 1 [H] , and upon thin¬ 
ning of the Riemann zero data set with ^ = 0.6, for all displayed values of s. This 
is in contrast to what was found in Figure [^ for the two-point function, where the 
accuracy deteriorates after approximately one and a half units of the mean spacing. 
A significant difference between the two quantities that may explain this observa¬ 
tion relates to the respective large s forms. That is, the spacing distribution decays 
as an exponential or faster (see (3.41) with (3.44)-(3.47)), whereas the correction 
term for the two-point correlation function is an oscillatory order one quantity for 
all s (see (2.2)). It follows that at a graphical level discrepancies in the case of 
the correction term to the spacing distribution functional formal form cannot be 
quantitatively probed in distinction to the situation with the correction term for 
the two-point correlation function. 

In conclusion, our study corroborates the study of [S], and so lends further 
weight to the conjecture that correction terms to the Montgomery-Odlyzko law 























20 


PETER J. FORRESTER AND ANTHONY MAYS 



Figure 9. Comparison of Riemann zero nearest neighbour spacing 
(where each zero is scaled by the leading term in O) with the 
corresponding DE solution {ie. the second derivative of ( |3.16 )). 
We also compare these to an extrapolated limit of a sequence of 
finite calculations as in (3.11). 




Figure 10. Comparison of next-to-leading order corrections to 
the quantities in Figure]^ 


themselves have a random matrix interpretation. Specifically, as put forward in 
[8] as an extension of the earlier work |39], the results of our study are consistent 
with the hypothesis that the correction terms to the Montgomery-Odlyzko law for 
correlation functions and associated distribution functions of the Riemann zeros 
coincides with the 0(1/1V^) correction terms for the corresponding quantities of the 
eigenvalues of random unitary matrices. This holds with the relationship between 
E and N given by (2.5), and with the change of scale (2.6). 
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